The universality of synchrony: critical behavior in a discrete model of stochastic 

phase coupled oscillators 
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We present the simplest discrete model to date that leads to synchronization of stochastic phase- 
coupled oscillators. In the mean field limit, the model exhibits a Hopf bifurcation and global 
oscillatory behavior as coupling crosses a critical value. When coupling between units is strictly 
local, the model undergoes a continuous phase transition which we characterize numerically using 
finite-size scaling analysis. In particular, the onset of global synchrony is marked by signatures of 
the XY universality class, including the appropriate classical exponents (3 and v, a lower critical 
dimension di c — 2, and an upper critical dimension d uc = 4. 
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In the early 1960's, experiments with the Belousov- 
Zhabotinsky reaction created a sensation by showing that 
dissipative structures and self-organization in systems far 
from equilibrium correspond to real observable physical 
phenomena. Since then, the breaking of time transla- 
tional symmetry has been a central theme in the analy- 
sis of nonlinear nonequilibrium systems. However, in the 
later studies of spatially distributed systems, most of the 
interest shifted to pattern forming instabilities, and sur- 
prisingly little attention was devoted to the question of 
bulk oscillation and the required spatial frequency and 
phase synchronization. On the other hand, the emer- 
gence of phase synchronization in populations of globally 
coupled phase oscillators, with the synchronous firing of 
fireflies as one of the spectacular examples, did gener- 
ate intense interest pj. Because intrinsically oscillating 
units with slightly different eigenfrequencies underlie the 
macroscopic behavior of an extensive range of biological, 
chemical, and physical systems, a great deal of litera- 
ture has focused on the mathematical principles govern- 
ing the competition between individual oscillatory ten- 
dencies and synchronous cooperation While most 
studies have focused on globally coupled units, leading 
to a mature understanding of the mean field behavior of 
several models, relatively little work has examined pop- 
ulations of oscillators in the locally coupled regime |jg. 
The description of emergent synchrony has largely been 
limited to small-scale and/or globally-coupled determin- 
istic systems Q, despite the fact that the dynamics of 
the physical systems in question likely reflect a combi- 
nation of finite-range forces and stochasticity. Two re- 
cent studies by Risler et al. Q represent notable excep- 
tions to this trend. They provide analytical evidence 
that locally-coupled identical noisy oscillators belong to 
the XY universality class, though to date there had been 
no empirical verification, numerical or otherwise, of their 
predictions. 

The difficulty with existing models of locally coupled 



oscillators is that each is typically described by a nonlin- 
ear differential equation, and it is notoriously computa- 
tionally intensive to deal with systems of coupled nonlin- 
ear differential equations, especially if they also involve 
a stochastic component. However, following Landau the- 
ory macroscopically observable changes occur with- 
out reference to microscopic specifics, instead giving rise 
to classes of universal behavior whose members may dif- 
fer greatly at the microscopic level. With this in mind, 
we construct the simplest model with short-ranged in- 
teractions between individual, stochastic, discrete phase 
units exhibiting global phase synchrony and amenable to 
extensive numerical study. 

Our starting point is a three-state unit |(j governed 
by transition rates g (Fig. Loosely speaking, we in- 
terpret the state designation as a generalized (discrete) 
phase, and the transitions between states, which we con- 
struct to be unidirectional, as a phase change and thus 
an oscillation of sorts. The probability of going from the 
current state i to state i + 1 in an infinitesimal time dt 
is gdt, with i— 1,2,3 modulo 3. For an isolated unit, the 
transition rate is simply a constant (g) that sets the os- 
cillator's intrinsic frequency; for many coupled units, we 
will allow the transition rate to depend on the neighbor- 
ing units in the spatial grid, thereby coupling neighboring 
phases. 




FIG. 1: Three state unit with generic transition rates g. 
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For an isolated unit we write the linear evolution equa- 
tion dP(i)/dt — MP(t), where the components P»(i) of 
the column vector P(t) = (P x (t) P 2 {t) P 3 (t)) T are the 
probabilities of being in state i at time t, and 



M = 




(1) 



The system reaches a steady state for Pj* = P 2 * = P3 = 
1/3. The transitions i — > i + 1 occur with a rough pe- 
riodicity determined by 3; that is, the time evolution of 
our simple model qualitatively resembles that of the dis- 
cretized phase of a generic noisy oscillator. 

The interesting behavior emerges when the transition 
probability of a given unit to the state ahead depends 
on the states of the unit's nearest neighbors in a spatial 
grid. To capture the physical nature of synchronization, 
we choose a function which compares the phase at a given 
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FIG. 2: (a) Onset of synchronization in d—3. The system 
size L=80 is used. Upper left inset: Fluctuations peak near 
the critical point, giving an estimation of a c =2.345±0.005. 
Pi and Pi undergo smooth temporal oscillations for large a 
(upper right inset), while lower a decreases temporal coher- 
ence (lower left inset). Lower right inset: Log-log plots of r 
vs L _1 with a=2.275, 2.3, 2.325, 2.375, 2.4, 2.425 from lowest 
to highest plots, (b) Finite size scaling analysis for d=3 us- 
ing the XY and Ising critical exponents. Data collapse with 
a c =2.345. 



site with its neighbors, and adjusts the phase at the given 
site so as to facilitate phase coherence. With universality 
in mind, we stress that the specific nature of the cou- 
pling is not important so long as we ultimately observe 
a transition to global synchrony at some finite value of 
the coupling parameter. For any unit, the transition rate 
from state i to state j is given by 



gij = gexp 



a{N 3 - Ni) 



2d 



(2) 



where the constant a is the coupling parameter and 5 is 
the Kronecker delta. is the number of nearest neigh- 
bors in state k, and 2d is the total number of nearest 
neighbors in d dimensional cubic lattices. While this 
choice is by no means unique and these rates are some- 
what distorted by their independence of the number of 
nearest neighbors in state i — 1, the form J2J is simpli- 
fied by this assumption and, as we shall see, does lead to 
synchronization. 

To test for the emergence of global synchrony, we first 
consider a mean field version of the model. In the large 
N limit with all-to-all coupling we write 



9ij = 5-exp [a(Pj - Pi)} S j>i+1 . 



(3) 



Note that in the mean field limit gtj does not depend on 
the location of the unit within the lattice. Also, there 
is an inherent assumption that we can replace Nk/N 
with Pfe. With this simplification we arrive at a nonlin- 
ear equation for the mean field probability, dP(t)/dt = 
M[P(t)]P(t), with 



-.912 .931 
M[P(t)} = \ .912 - 323 

323 -331, 



(4) 



Normalization allows us to eliminate P3 and obtain 
a closed set of equations for Pj and P%. We can further 
characterize the mean field solutions by linearizing about 
the fixed point (P*,P 2 *) = (1/3, 1/3). The complex con- 
jugate eigenvalues of the Jacobian evaluated at the fixed 
point, A± = g(2a — 3 ± zv / 3)/2, cross the imaginary axis 
at a = 1.5, indicative of a Hopf bifurcation at this value, 
which following a more detailed analysis can be shown 
to be supercritical. Hence, as a increases, the mean field 
undergoes a qualitative change from disorder to global 
oscillations, and the desired global synchrony emerges. 
Numerical solutions confirm this behavior, yielding re- 
sults that agree with simulations of an all-to-all coupling 
array ||. Here we characterize the breakdown of the 
mean field description for the nearest-neighbor coupling 
model as spatial dimension is decreased. 

We perform simulations of the locally coupled model 
in continuous time on d-dimensional cubic lattices with 
periodic boundary conditions. Time steps are 10 to 100 
times smaller than the fastest local average transition 
rate, i.e., dt <§C e~ a (we set 3=1). We find that much 
smaller time steps lead to essentially the same results. 
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Starting from random initial conditions, all simulations 
were run until an apparent steady state was reached, and 
statistics are based on 100 independent trials. 

To characterize the emergence of phase synchrony, we 
introduce the order parameter 



(R), R 



1 

N 



N 



(5) 



Here cf>j is the discrete phase 27r(fc-l)/3 for state k € 
{1,2,3} at site j. The brackets represent an average 
over time in the steady state and over all independent 
trials. Nonzero r in the thermodynamic limit indicates 
synchrony. We also calculate the generalized susceptibil- 
ityx-^Ki? 2 )-^) 2 ]. 

In d=2 we do not see the emergence of global oscilla- 
tory behavior. Instead, we observe intermittent oscilla- 
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FIG. 3: Transition in d—4 (top) and d=5 (bottom): The or- 
der parameter near the transition point is shown for various 
system sizes. O L=4, □ L=8, x L=12, A L=16 for d=4 and 
O L=i, □ L=6, x L=8, A £=10 for d=5. The upper left 
inset in each panel shows the generalized susceptibility which 
peaks at a=1.900±0.025 for d=A and at a=1.750±0.015 for 
d—5. The lower right inset shows the system size dependence 
of the order parameter. For d=4 the coupling constant varies 
from a=1.6 to 2.4 in increments of 0.1 (excluding o=1.9) from 
lowest to highest plots and for d=5 from a— 1.4 to 2.2 in in- 
crements of 0.1. 



tions (for very large values of a) that decrease drastically 
with increasing system size. In fact, r — * in the ther- 
modynamic limit, even for very large values of a |||. We 
conclude that the phase transition to synchrony cannot 
occur for d—2. Interestingly, snapshots of the system re- 
veal increased spatial clustering as a is increased, as well 
as the presence of defect structures, perhaps indicative of 
Kosterlitz-Thouless-type phenomena 5] . Further studies 
along these lines are underway. 

In contrast to the d=2 case, which serves as the lower 
critical dimension, a clear thermodynamic-like phase 
transition occurs in three dimensions. We see the emer- 
gence of global oscillatory behavior as a increases past 
a critical value a c . Figure 0a, shows the behavior of the 
order parameter as a is increased for the largest system 
studied (L=80); the upper left inset shows the peak in x 
at a=2.345±0.005, thus providing an estimate of the crit- 
ical point a c . We see no change as system size is increased 
beyond L=40. At any rate, finite size effects are within 
the range of our estimation. The lower right inset in 
Fig-Et shows explicitly that for a < a c , r — * as system 
size is increased, and a disordered phase persists in the 
thermodynamic limit. For a > a Cl the order parameter 
approaches a finite value as the system size increases. We 
tried to apply the Binder cumulant crossing method 
for determining a c more precisely, but residual finite size 
effects and statistical uncertainties in the data prevent 
us from determining the crossing point with more preci- 
sion than that stated above. In any case, the accuracy 
of our current estimation of the critical point suffices to 
determine the universality class of the transition. 

To further characterize this transition, we use finite 
size scaling analysis by assuming the standard scaling 



= L~F[(a - a c )L-]. 



(6) 



Here F(x) is a scaling function that approaches a con- 
stant as x — ► 0. To test our numerical data against differ- 
ent universality classes we choose the appropriate critical 
exponents for each, recognizing that there are variations 
in the reported values of these exponents 0. For the 
XY universality class we use the exponents /3=0.34 and 
j/=0.66 [lJ. For the Ising exponents we use /3=0.31 and 
^=0.64 [l^. In Fig. [2b, we see quite convincingly a col- 
lapse when exponents from the XY class are used. For 
comparison, we also show the data collapse with 3D Ising 
exponents (note the scale differences). 

For g?=4 we estimate the transition coupling to be 
a c =1.900±0.025 from the peak in \ (Fig-Efc)- Because 
we expect (i=4 to be the upper critical dimension in ac- 
cordance with XY/Ising behavior, we anticipate a slight 
breakdown of the scaling relation A priori it is not 
clear how strongly JBJ will be violated in d=4. As shown 
in Fig. the data collapse is very good with the mean 
field exponents. As such, our simulations suggest that 
d=4 serves as the upper critical dimension; additionally, 
it appears that corrections to finite-size scaling at d—4 
are not substantial, though a much more precise study 



4 



o L-4 

a L-8 

« L=l 2 

A L-16 



20 -10 10 20 
(a/a - 1)L ,A ' 



A 4. 



-30 



-20 



10 10 

(a/a c -l)L ,/v 



O L=4 

□ L=6 

* L=8 

A L=10 

20 



5 

4.5 
4 

3.5 
3 4 
2.5 • 
2 

1.5 
1 

0.5 




FIG. 4: Finite size scaling analysis for d=4 and d=5: Data 
collapse using ansatz JfoJ with mean field exponents. 



would be needed to investigate such corrections in greater 
detail. 

To further support the claim that d uc =4, we con- 
sider the case d=5. We see a transition to synchrony 
at a c =1.750±0.015 (Fig. |3fc>). As expected, this value 
is considerably closer than the critical coupling in four 
dimensions to the mean field value a c =1.5. The data 
collapse with the mean field exponents is excellent, as 
shown in Fig. We note the rarity of computations in 
such a high dimension. 

In conclusion, while nonequilibrium phase transitions 
exhibit a much wider diversity in universality classes than 
equilibrium ones 01 > ^ is remarkable that the prototype 
of a nonequilibrium transition, namely, a phase transi- 
tion that breaks the symmetry of translation in time, is 



described by an equilibrium universality class. In par- 
ticular, the Mermin- Wagner theorem, stating that con- 
tinuous symmetries can not be broken in dimension two 
or lower, appears to apply. The XY model is known 
to display a Kosterlitz-Thouless transition in which, be- 
yond a critical temperature, vortex pairs can unbind into 
individual units creating long range correlations. Prelim- 
inary results indicate that a similar transition occurs in 
our model. 

Finally, a note of caution concerning the discreteness 
of the phase is in order. We first note that microscopic 
models often feature discrete degrees of freedom. For ex- 
ample, our model is reminiscent of the triangular reaction 
model of Onsager 0] , on the basis of which he illustrated 
the concept of detailed balance as a characterization of 
equilibrium. Continuous phase models appear in a suit- 
able thermodynamic limit. We stress that the breaking 
of time translational symmetry can occur independently 
of whether the phase is a discrete or continuous variable. 
It is, however, not evident whether continuous and dis- 
crete phase models belong to the same universality class. 
For example, the three state ferromagnetic Potts model 
displays a weak first order phase transition in d=3 
while the anti-ferromagnetic version belongs to the XY 
universality class |lCt [16( . The results found here appear 
to be compatible with the latter, but a renormalization 
calculation confirming this hypothesis would be welcome. 
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